# =============================================================================
# BOUNDARY CONDITION PERFORMANCE BENCHMARKS
# =============================================================================

using BenchmarkTools
using PDEJulia
using PDEJulia.Grids
using LinearAlgebra
using Printf

"""
Comprehensive benchmarking suite for boundary condition performance.
Tests various BC types, grid sizes, and application scenarios.
"""

function run_boundary_benchmarks()
    println("=" ^ 80)
    println("BOUNDARY CONDITION PERFORMANCE BENCHMARKS")
    println("=" ^ 80)
    println()
    
    # Test grid sizes
    small_size = 100
    medium_size = 1000  
    large_size = 10000
    
    # Test different BC types
    bc_types = [
        ("Dirichlet", DirichletBC(1.0)),
        ("Neumann", NeumannBC(0.5)), 
        ("Robin", RobinBC(1.0, 0.1, 2.0)),
        ("Periodic", PeriodicBC()),
        ("Mixed", MixedBC(DirichletBC(1.0), NeumannBC(0.0))),
        ("Time-Dependent", TimeDependentBC(t -> DirichletBC(sin(t))))
    ]
    
    println("1. BC APPLICATION PERFORMANCE")
    println("-" ^ 50)
    
    for (name, bc) in bc_types
        println("\\n$name Boundary Conditions:")
        benchmark_bc_application(bc, small_size, medium_size, large_size)
    end
    
    println("\\n\\n2. OPERATOR PERFORMANCE WITH BCs") 
    println("-" ^ 50)
    benchmark_operator_performance()
    
    println("\\n\\n3. MEMORY ALLOCATION ANALYSIS")
    println("-" ^ 50)
    benchmark_memory_allocation()
    
    println("\\n\\n4. 2D/3D SCALING TESTS")
    println("-" ^ 50)
    benchmark_multidimensional()
    
    println("\\n\\n5. COMPARISON WITH py-pde EQUIVALENT")
    println("-" ^ 50)
    benchmark_vs_py_pde_patterns()
    
    println("\\n" * "=" ^ 80)
    println("BENCHMARK SUMMARY COMPLETE")
    println("=" ^ 80)
end

"""
Benchmark boundary condition application for different grid sizes.
"""
function benchmark_bc_application(bc::BoundaryCondition, small_size::Int, medium_size::Int, large_size::Int)
    # Test 1D grids
    for size in [small_size, medium_size, large_size]
        grid = UnitGrid(size, (0.0, 1.0))
        field_data = rand(size)
        
        set_boundary_conditions!(grid, [bc])
        
        # Benchmark application
        result = @benchmark apply_boundary_conditions!($field_data, $grid) setup=(field_data = rand($size))
        
        avg_time = mean(result.times) / 1e6  # Convert to ms
        allocs = result.allocs
        memory = result.memory
        
        @printf "  Size %6d: %8.3f ms  (%d allocs, %s)\\n" size avg_time allocs Base.format_bytes(memory)
    end
end

"""
Benchmark differential operators with boundary conditions.
"""
function benchmark_operator_performance()
    grid_1d = UnitGrid(1000, (0.0, 1.0))
    grid_2d = CartesianGrid([100, 100], [(0.0, 1.0), (0.0, 1.0)])
    
    field_1d = rand(1000)
    field_2d = rand(100, 100)
    
    bc_types = [
        ("Dirichlet", DirichletBC(0.0)),
        ("Neumann", NeumannBC(0.0)),
        ("Periodic", PeriodicBC())
    ]
    
    println("\\n1D Laplacian Performance:")
    for (name, bc) in bc_types
        set_boundary_conditions!(grid_1d, [bc])
        
        result = @benchmark adaptive_laplace($grid_1d, $field_1d) setup=(field_1d = rand(1000))
        avg_time = mean(result.times) / 1e6
        @printf "  %-12s: %8.3f ms\\n" name avg_time
    end
    
    println("\\n2D Laplacian Performance:")
    for (name, bc) in bc_types
        set_boundary_conditions!(grid_2d, [bc, bc])
        
        result = @benchmark adaptive_laplace($grid_2d, $field_2d) setup=(field_2d = rand(100, 100))
        avg_time = mean(result.times) / 1e6
        @printf "  %-12s: %8.3f ms\\n" name avg_time
    end
end

"""
Analyze memory allocation patterns for different BC operations.
"""
function benchmark_memory_allocation()
    grid = UnitGrid(1000, (0.0, 1.0))
    field_data = rand(1000)
    
    bc_types = [
        ("Dirichlet", DirichletBC(1.0)),
        ("Neumann", NeumannBC(0.5)),
        ("Robin", RobinBC(1.0, 0.1)),
        ("Periodic", PeriodicBC())
    ]
    
    println("\\nMemory Allocation Analysis (1000 points):")
    for (name, bc) in bc_types
        set_boundary_conditions!(grid, [bc])
        
        # Measure allocations
        field_copy = copy(field_data)
        allocs_before = Base.gc_alloc_count()
        apply_boundary_conditions!(field_copy, grid)
        allocs_after = Base.gc_alloc_count()
        
        # Benchmark for detailed memory info
        result = @benchmark apply_boundary_conditions!($field_copy, $grid) setup=(field_copy = copy($field_data))
        
        @printf "  %-12s: %d allocs, %s\\n" name result.allocs Base.format_bytes(result.memory)
    end
end

"""
Test scaling behavior with multidimensional grids.
"""
function benchmark_multidimensional()
    # 2D grids of increasing size
    sizes_2d = [32, 64, 128]
    
    println("\\n2D Grid Scaling (Dirichlet BC):")
    for size in sizes_2d
        grid = CartesianGrid([size, size], [(0.0, 1.0), (0.0, 1.0)])
        field_data = rand(size, size)
        
        bcs = [DirichletBC(0.0), DirichletBC(0.0)]
        set_boundary_conditions!(grid, bcs)
        
        result = @benchmark apply_boundary_conditions!($field_data, $grid) setup=(field_data = rand($size, $size))
        avg_time = mean(result.times) / 1e6
        
        @printf "  %dx%d grid: %8.3f ms\\n" size size avg_time
    end
    
    # 3D test (smaller sizes due to memory)
    if false  # Enable for full 3D testing
        println("\\n3D Grid Test:")
        grid_3d = CartesianGrid([32, 32, 32], [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)])
        field_3d = rand(32, 32, 32)
        
        bcs_3d = [DirichletBC(0.0), DirichletBC(0.0), DirichletBC(0.0)]
        set_boundary_conditions!(grid_3d, bcs_3d)
        
        result = @benchmark apply_boundary_conditions!($field_3d, $grid_3d) setup=(field_3d = rand(32, 32, 32))
        avg_time = mean(result.times) / 1e6
        
        @printf "  32x32x32 grid: %8.3f ms\\n" avg_time
    end
end

"""
Compare with patterns typical in py-pde workflows.
"""
function benchmark_vs_py_pde_patterns()
    # Simulate typical py-pde usage patterns
    
    println("\\nTypical PDE Workflow Patterns:")
    
    # Pattern 1: Diffusion with Dirichlet BCs
    grid = UnitGrid(256, (0.0, 1.0))
    field = rand(256)
    bcs = [DirichletBC(0.0)]
    set_boundary_conditions!(grid, bcs)
    
    # Time a full diffusion step (BC application + Laplacian)
    result = @benchmark begin
        apply_boundary_conditions!($field, $grid)
        adaptive_laplace($grid, $field)
    end setup=(field = rand(256))
    
    avg_time = mean(result.times) / 1e6
    @printf "  Diffusion step (1D, 256 pts): %8.3f ms\\n" avg_time
    
    # Pattern 2: 2D reaction-diffusion with mixed BCs
    grid_2d = CartesianGrid([64, 64], [(0.0, 1.0), (0.0, 1.0)])
    field_2d = rand(64, 64)
    bcs_2d = [PeriodicBC(), NeumannBC(0.0)]
    set_boundary_conditions!(grid_2d, bcs_2d)
    
    result_2d = @benchmark begin
        apply_boundary_conditions!($field_2d, $grid_2d) 
        adaptive_laplace($grid_2d, $field_2d)
    end setup=(field_2d = rand(64, 64))
    
    avg_time_2d = mean(result_2d.times) / 1e6
    @printf "  RD step (2D, 64x64): %8.3f ms\\n" avg_time_2d
    
    # Pattern 3: Time-dependent boundary
    td_grid = UnitGrid(128, (0.0, 1.0))
    td_field = rand(128)
    td_bc = TimeDependentBC(t -> DirichletBC(sin(t)))
    set_boundary_conditions!(td_grid, [td_bc])
    
    result_td = @benchmark apply_boundary_conditions!($td_field, $td_grid; time=1.5) setup=(td_field = rand(128))
    avg_time_td = mean(result_td.times) / 1e6
    @printf "  Time-dependent BC: %8.3f ms\\n" avg_time_td
end

"""
Run performance profiling to identify bottlenecks.
"""
function profile_boundary_performance()
    println("\\nPROFILING BOUNDARY CONDITION PERFORMANCE")
    println("-" ^ 50)
    
    using Profile
    
    # Set up test case
    grid = CartesianGrid([100, 100], [(0.0, 1.0), (0.0, 1.0)])
    field_data = rand(100, 100)
    bcs = [MixedBC(DirichletBC(1.0), NeumannBC(0.0)), RobinBC(1.0, 0.1, 2.0)]
    set_boundary_conditions!(grid, bcs)
    
    # Profile the operation
    Profile.clear()
    @profile for i in 1:1000
        field_copy = copy(field_data)
        apply_boundary_conditions!(field_copy, grid)
        adaptive_laplace(grid, field_copy)
    end
    
    # Show profile results
    Profile.print(mincount=10)
end

"""
Test boundary condition accuracy vs performance tradeoffs.
"""
function benchmark_accuracy_vs_performance()
    println("\\nACCURACY vs PERFORMANCE ANALYSIS")
    println("-" ^ 50)
    
    # Test different finite difference orders
    grid = UnitGrid(101, (0.0, π))
    
    # Analytical test case: sin(x) with known Laplacian -sin(x)
    x = [π * i / 100 for i in 0:100]
    analytical_field = sin.(x)
    analytical_laplacian = -sin.(x)
    
    # Test with Dirichlet BCs that match analytical solution
    bcs = [DirichletBC(0.0)]  # sin(0) = sin(π) = 0
    set_boundary_conditions!(grid, bcs)
    
    # Benchmark and measure accuracy
    result = @benchmark adaptive_laplace($grid, $analytical_field)
    computed_laplacian = adaptive_laplace(grid, analytical_field)
    
    error = norm(computed_laplacian - analytical_laplacian, Inf)
    avg_time = mean(result.times) / 1e6
    
    @printf "  Laplacian accuracy test:\\n"
    @printf "    Max error: %12.2e\\n" error
    @printf "    Time:      %12.3f ms\\n" avg_time
    @printf "    Efficiency: %11.2e error/ms\\n" error/avg_time
end

# Run all benchmarks when script is executed
if abspath(PROGRAM_FILE) == @__FILE__
    run_boundary_benchmarks()
    # profile_boundary_performance()  # Uncomment for profiling
    benchmark_accuracy_vs_performance()
end